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Abstract 

A fast method of an arbitrary high order for approximating volume potentials is 
proposed, which is effective also in high dimensional cases. Basis functions intro- 
duced in the theory of approximate approximations are used. Results of numerical 
experiments, which show approximation order 0(h a ) for the Newton potential in high 
dimensions, for example, for n = 200 000, are provided. The computation time scales 
linearly in the space dimension. New one-dimensional integral representations with 
separable integrands of the potentials of advection-diffusion and heat equations are 
obtained. 



1 Introduction 



The cubature of high-dimensional volume potentials plays an important role in a wide 
range of applications in physics, chemistry, biology, financial mathematics etc. Even a few 
years ago this problem encountered unsurmountable difficulties due to the so-called "curse 
of dimensionality" . With the development of separated representations (also called tensor- 
structured approximations) by Beylkin and Mohlenkamp in O H] the problem became 
tractable. In fact, in recent years several fast algorithms for the computation of multi- 
dimensional convolutions with singular kernels have been proposed, cf. [21 [61 El [5] . 
They are mainly based on piecewise polynomial approximations of a separated repre- 
sentation of the density. Then a suitable separated approximation of the action of the 
convolution operator on the basis functions allows one to determine the multi-dimensional 



1 



convolution in question by computing a number of one-dimensional convolutions. On this 
way, the complexity of computing the multi-dimensional convolutions on a uniform tensor- 
product grid of size h can be reduced from 0(h~ n \ \ogh\), achieved with the traditional 
multi-variate FFT, to 0{nh~ 1 \ log h\), where n is the space dimension. A certain drawback 
of the schemes used in present is the necessity to find accurate separated representations 
of the convolution operator acting on piecewise polynomials. This procedure has been 
developed only for a few particular kernels ([21 Ej) and is rather involved, especially for 
higher order approximations. 

In this article we propose a fast method of an arbitrary high order for approximating 
volume potentials, which is effective also in high dimensional cases. We use basis func- 
tions introduced in the theory of approximate approximations I12j. see also |13j and 
the references therein. We report on numerical experiments which show approximation 
order 0(h s ) for the Newton potential up to dimension n = 200 000. The computational 
complexity of the algorithm scales linearly in the physical dimension. In the experiments, 
high order cubatures lead to a minor increase of the computation time in comparison with 
second order approximations, but are some orders of magnitude more accurate, especially 
in high dimension. 

The method combines our basis functions with an approach from Khoromskij [9] for 
computing volume potentials on uniform or composite refined grids. For brevity of ex- 
position we restrict ourselves to the case of uniform grids with step size h. In [9], the 
density is approximated by piecewise constants and the action of the convolution operator 
on the basis function is written as a one-dimensional integral with a separable integrand, 
i.e., a product of functions depending only on one of the space variables. Then an accu- 
rate quadrature rule of this one-dimensional integral provides a separated representation 
of the convolution operator. The convergence orders 0{h 2 ) and 0(/i 3 ) using Richardson 
extrapolation are confirmed by numerical experiments for n = 3. 

Instead of piecewise constants we use basis functions, which are Gaussians or products 
of Gaussians and special polynomials and give rise to high-order semi-analytic cubature 
formulas for volume potentials. In Section [2] we describe these formulas and give error 
estimates. According to [H] (see also [131 Section 6.3]), the action of volume potentials 
on the basis functions allows for one-dimensional integral representations with separable 
integrands. In Section [3] we demonstrate this by the example of Newton's potential. 
We also obtain new one-dimensional representations for the potentials of the advection- 
diffusion and the heat equations. These one-dimensional integrals in combination with a 
quadrature rule lead to accurate separated representations of the potentials. In Sectional 
we provide results of numerical experiments for the Newton potential showing that even for 
very high space dimensions these approximations preserve the predicted convergence order 
of the cubature. In the final section [5] we describe the quadrature of the one-dimensional 
integral representations used in the numerical examples. Here we also report on tests for 
second and fourth order formulas of the Newton potential and the inverse of — A+a 2 , which 
provide estimates of the rank of the separated representation required to approximate the 
action of the potential on a basis function with a prescribed relative error. 

2 Semi-analytic cubature formulas for potentials 

Here we collect some results on high-order cubature formulas for the volume potentials of 
the differential operators —A and —A + 2b • V + c and of the heat potential in W 1 . More 
details can be found in [13]. 
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2.1 Newton potentials in IR n 

The Newton potential is the inverse of the Laplace operator and given in M. n by the formula 

C nU ( x ) = ~ l) [ " (y) - dy, n > 3 . (2.1) 
n y 1 4W 2 J |x-y| n " 2 

Here and in the following we denote by bold x a vector of length n with the components 
x = (x±, ... ,x n ). The integral (|2.ip is a unique solution of Poisson's equation 

-Af = u in IT, |/(x)| < C\x\ n ~ 2 as |x| -> oo . 

In [llj IS] cubature formulas for computing (|2.ip were constructed which are based on the 
approximation of the density u by functions with analytically known Newton potentials. 
In particular, one can approximate u by the quasi-interpolant 

ti A (x) = P-"/ 2 £ u(/im) r? 2M » (2-2) 



with the basis function 



r?2M (x) = 7r-"/ 2 4y_ 2 1 ) (|x| 2 )e-l x l 2 , MEN, (2.3) 



where are the generalized Laguerre polynomials 



4-"«^a)'(e-»,-), 7>-l- 



The Newton potential of r\2M is given by 



WW) " iW^ 7 ^ " »• W) + V ? '4 (J + 1) (2 ' 4) 



J=0 

with the incomplete Gamma function 



7M = /^.-*. 





Then the sum 



ft- 2 \ - /x — /im\ 

£nUh{x) = pw /2-i 2^ nhm)L n rj2M [ ) (2.5) 



m£Z 

is a semi-analytic cubature formula for It has been shown, that for sufficiently smooth 
and compactly supported functions (|2.5|) approximates C n u with the error 

0(h 2M ) + 0(e"^ 2 ft 2 ) . (2.6) 

This follows from the general asymptotic error estimate 0(h 2M )+e for the quasi-interpolant 
(|2.2p with sufficiently smooth and decaying generating functions r\ satisfying the moment 
condition 

J ??(x) dx = 1 , J x a r](x) dx = 0, V a , 1 < \a\ < 2M , (2.7) 
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see [13\ Chapter 2]. The saturation error e does not converge to zero as h —* 0, however, 
because of 



e = 0[max\ ^ Trj(VVk)e Zm ^' x)/n ), (2. 
keZ"\{o} 



x 



it can be made arbitrarily small if the parameter T> is sufficiently large. Additionally, the 
Newton potential maps the fast oscillating saturation term (|2.8p into a function with norm 
of order 0(h 2 s), which establishes the estimate ()2.6p . Therefore, in numerical computa- 
tions with T> > 3 formula (|2.5|) behaves like a high order cubature formula. 

2.2 Potentials of advection-diffusion operators in M n 

The fundamental solution of the operator — A + 2b • V + c with a vector b £ C n and c E C 
depends on the value of c + |b| 2 . Here we use the notation 

n 

(y,z) = ^VjZj and |y| 2 = (y,y) 
i=i 

also for vectors y, z G C™. If c + |b| 2 7^ 0, then the fundamental solution can be given as 

e ( b > x ) /IxK 1 -™/ 2 

KAW= (2^7HtJ ^n/2-l(A|x|), 

where A 6 C \ (- 00, 0] with A 2 = c + |b| 2 and K v is the modified Bessel function of the 
second kind, also known as Macdonald function, [1, 9.6]. If c + |b| 2 = 0, then for n > 3 

«o(x) 



T(§ - 1) e < b ' x > 



47T n / 2 |x| n ~ 2 

To derive a cubature formula for the volume potential with the kernel k,\ we look for a 
solution of 

- A/ + 2b- V/ + c/ = e" |x|2 , x£ R n , (2.9) 
which is given as the one-dimensional integral 



-|x| 2 A n/2-l /• 

/(x) = |2x + b|n/2 _ 1 J ^n/2- 1 (Ar)/ n/2 _ 1 (2|x + b|r)r e^ dr, (2.10) 


where I y is the modified Bessel function of the first kind, see \13\ Section 5.2]. 

Using known analytic expressions of I n+ \/2 and K n+ i/ 2 (cf. [I]) it is possible to derive 
analytic formulas of (|2.10p for odd space dimensions n. In particular, if n = 3, then 

7r e~' x ' 2 / / i , . , . ,>\ (i 



/W = X^TM I" (2 < A - I2x + b|) j - w (-(A + |2x + b| 
where w denotes the scaled complementary error function 

z 

w ( z ) = e ~ z2 erfc(-iz) = e" 22 (l + j e' 2 dtj 
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and erfc = 1 — erf is the complementary error function. From 

«(W 2 )e- |xp =E^A fc e-N 2 (2.11) 

A:=0 

(see [131 Theorem 3.5]), one can derive as in the case of Newton potentials semi-analytic 
cubature formulas for the potential of the advection-diffusion equation with the approxi- 
mation rate (|2,6p . 



2.3 Heat potential 

Next we consider the non-homogeneous heat equation 

ft - z^A x / = «(x, t) t > 0, v > , 



Dll 



(2.12) 



/(x,0)=0 xG 

It is well known that the solution of this Cauchy problem can be written as 

t 

/(x,t) = yVt-A<,A))(x) dX, 
o 

where P| is the Poisson integral 

(Pt«(-,A))(x) = / e-l— l 2 /(-) u(y , A)dy . 

Suppose that the right-hand side u is a sufficiently smooth function with compact support 
in W 1 x An approximation of the solution f(x,t) of (|2.12p can be obtained if u is 
approximated by the quasi-interpolant on the rectangular grid {(hm, ri)} 

UhAX} t ) = n ~^—^ Yl < hm ' ™) *- (t - Tl)2/{VoT2) e-\*-W/<PW . (2.13) 

yJV G V n 

It can be easily seen, that 

h n e -(\-ri) 2 /(V T 2 ) e -\x-h. m \ 2 /(Vh 2 +4u(t-\)) 
(P t -A^(-,A))(x) = -^= £ «(^") (2>tf + M t_ A ))«/2 ' 

meZ" 

hence the sum 

t 

(7 ? t -AMfe,r(-,A))(x)dA = — y— == V K 2 (x- hm,t,Ti)u(hm,Ti) 



o 

u mez™ 
is an approximation of /(x, t), where we use the notation 

} (t-\-Tif/(T>oT 2 ) -\ X \*/(ph 2 +4v\) 
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Since it^ T (x, t) differs from u(x,t) by 

|ti(x, t) - «/, , T (x, t)| < e + c((r V / Ai) 2 + (/iv 7 ^) 2 ) , Vx G M n , t G [0, T] , 
the function //i )T (x, t) approximates the solution /(x, i) with the error 

i r /■ e -l x -yl 2 /(4(*-A)) 

|/(x,t)-/ hiT (x,t)| = — -— 2 y rfA y (t _ A)w/2 Ky,A)-^, r (y,A)|dy 



(4vr^) n / 2 7 y (t - A)"/ 2 

R» 

<T\\u- u/ l) rlU oo (K n x[0,ri) < e + c((i"\/Po) 2 + {hVv) 2 ) 
for all x G M n and i G [0,T]. 

3 One-dimensional integral representations of the potentials 
acting on Gaussians 

The computation of the mentioned cubature formulas on the grid {hk} leads to discrete 
convolutions 

^2 a k-m u(hm) , (3.1) 

m 

where the indices m, k belong to some subset of Z n , and the coefficients are given 
either analytically or by a one-dimensional integral with smooth integrand. For general 
functions u the most efficient summation methods for (|3.ip are probably fast convolutions 
based on multi-variate FFT. However, even for the space dimension n = 3 problems of 
moderate size often exceed the capacity of available computer systems. 

The situation in high dimensions is much better when the vectors {u(hm)} and {a^} 
allow separated representations [3], i.e. for given accuracy e they can be represented as a 
sum of products of vectors in dimension 1 

R n 

u(hmi, . . . ,hm n ) = ^ r p |~| uf \hrrij) + 0(e) , 
p=i j=i 

R. n 

a(h, . . . , k n ) = J^sp II v f( k l) + °( £ ) • 
p=i j=i 



Then an approximate value of (|3.ip can be computed by the sum of products of one- 
dimensional convolutions 

R n 

^2 ° k - m u ( hm ) ~ r v s i II 5Z v i ( k 3 ~ m i) u T (hrrij) . 

m P,9=l i=l m j 

To obtain a separated representation of the vector {a^} the following idea from Khoromskij 
[S] can be applied: Suppose that the density u is approximated by interpolation or quasi- 
interpolation using linear combinations of the translates </>(• — hm) of a basis function. 
Then the components of {a^} are the values of the volume potential with the kernel g 
acting on <j>, 



a-k = J ff(k - y)<f>(y) dy 
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If this integral can be transformed to a one-dimensional integral with separable integrand 



°? n 

g(x - y)0(y) dy = J Y[9j(xj,t)dt 



3=1 

then a separated approximation of {a^} is obtained by applying an accurate quadrature 
rule for this integral 

R n 

«k ~ ue TT 9j(kj,Te) . 

1=1 3=1 

In [9] one-dimensional integral representations were obtained for the potentials of the 
Laplace and of the advection-diffusion equation with b = 0, c > 0, applied to the charac- 
teristic functions <j) of a cube in M 3 . This gives rise to a fast algorithm which is based on 
piecewise constant basis functions. The convergence order 0(h 2 ) is proved together with 
the order 0(h 3 ) by using Richardson extrapolation. 

In this section we describe one-dimensional integral representations with separable in- 
tegrands for the Newton potential and the potentials of the advection-diffusion and the 
heat equations acting on basis functions introduced in the previous section. 

3.1 Newton potential 

The separated representation of the second order cubature formula 

£ u(Am)4>(e-H a )(r m ) (3.2) 

for the Newton potential with 

x — hm , _ , _i.i2 w N 1 (n , , 2 



— 11111 / _i i2 \ . . ± / /t 

— = — and £ n e 1 1 (x) = — — : — ^7 1, x 

y/Vh V ' 4|x|™-2 'V2 1 

follows from the formula 

1 °f e -|x| 2 /(l+t) 



£ n (e~ M )(x) = - A I ° ' " ' , ' dt , 

V 4 7 (1 +t)"/2 



which was obtained in [H] and is valid for n > 3. The quadrature rule of the last integral 
w-|. Pv , if e-lxlVd+rp) i« JL e-^/d+r.) 

Me )M« 3 I> (1 + = 3 ^n7T+^)W 

p=i v p=i j=i v F/ 

with certain quadrature weights u; p and nodes r p gives the separated approximation 

m 2 * » e -fe-m,) 2 /^(l+^)) 

V ' m£Z" k=l j=l V fey 

To get an approximation for higher order cubature formulas we note that the same 
convergence order 0(/i 2M ) + 0(e-^ h 2 ) as in the case of generating functions (|2.3p holds, 
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when the density is approximated by the sum 

Mm<x) = (vrP)-/ 2 ]T u(hm) f[rj 2M fXj " ^ 

Here the basis function is the tensor product of the one-dimensional functions 

W*)=41 2 i(* 2 )e-* 2 , (3-3) 
which obviously satisfies the moment condition (|2.7p . Thus the cubature formula 

^ u{hra)C n ( JJ ~ 
meZ" i=i 

approximates C n u with the error (j2.6|) . The relation 



,2 n J, 

= o^=r E "H4(nHra (3-4) 



»(■) = E e E 4- 1/2 V), (3.5) 



leads to the one-dimensional integral representation of C n ( IIj=i V2M^j by writing the 
solution of the Poisson equation 



-Ati(x) = JJr/2M( 



as the integral 



1 JL ? e -|x| 2 /(W) ^ 

= i 7 ( TT V eeII! _e__ e -?/( 1 + i ) ^ dt 

4 y \11 ^ A;!4 fc dx 2fc /(1 + iW 2 

n .7=1 A:=0 J y 



o i =1 /«=o 

nM-1 _*»/(l+t) , ,,„ w T 2 



_1 f/^^ e^ ( - 1/2) / g xx g 

4 yviiz^ (1+*)* k \i + tJ) ri + tW2" 

Q j = l fc = ^ ' 

Hence the separated representation of the Newton potential is reduced to find accurate 
quadrature rules for the parameter dependent integrals 

00 n M—l 2 

*-w - / E Ir ^4-" 2 >( 1 ^)) mm. (3.6) 

It is necessary that the quadrature rule approximates the integrals (|3.6p with prescribed 
error for the parameters xj = (kj — m^/y/V within the range \xA < K and some given 
bound K. This will be discussed in Section [SJ 
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3.2 Potentials of advection-diffusion operators 

Here we look for a one-dimensional integral representation with separable integrand of 
the integral (|2.10p . Let / be the solution of the advection-diffusion equation (j2.9j) . Then 
v = /e~^ b ' x ^ satisfies 

—Av + X 2 v = e -l x+b / 2 ! 2 el b l 2 / 4 with A 2 = c + |b| 2 

in IR n . In was shown in |14j that v can be represented in the form 

e |b| 2 /4 °f e -A 2 i/4 -|x+b/2|V(l+t) 

«(x) = ^— — / j- -r dt (3.7) 

o 

for Re A 2 > and n > 3. If Re A 2 > 0, then (J37f]) is valid for all space dimensions n, see 
also |13} Theorem 6.4]. Thus 

r e (b,x) e |b| 2 /4 7 e -*( c+ l b l 2 )/ 4 e-l x+b / 2 l 2 /( 1+ *) 
/(x) = y« A (x-y)e-W dy = J dt 



oo 

1 r e~ ct / 2 

2 7(1 + 2t) n / 2 



|x-ib| 2 /(l+2t) dt (3 _ 8) 



Consequently, if Re(c+ |b| 2 ) > 0, then an approximate solution of the advection-diffusion 
equation 

-A/ + 2b • V/ + cf = u 

is given in M n by the sum 

Vh 2 ^ , 7 e" m2c */ 2 



A(x) = Vh y u{hm) f e C 7 e-lx-fan-^PhWCPfc'd+aO) dt . 
J v y 2(ttPW 2 ^ 7 (l + 2tW 2 

x ' ma" q x ' 



which approximates / with the order 0(h 2 ). Analogously to the case of Newton potentials 
the one-dimensional integrals for cubature formulas of order (|2.6p are based on the basis 
functions Y\ V2M (%j ) • Now the relation (|3.5j) leads to 



/ " ^ ( — 1) " f 2 

«A(x " y) II ^m(%) ^ = IIE IfF ^2fc / k a(x - y) e~l y l dy 

J — 1 J — 1 fc=0 ^ 

= - / e fTT V d -ix-«biVfi+2tA df 

27 (l + 2tW2Ul fc!4 fc dx 2fc e J ar 

oo 

= 2 / (7^p (lI»(2^i-%))^ (3-9) 

n .7=1 



with the function 



fc=0 v ; 
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3.3 Heat potential 

Recall that the approximation of the non-homogeneous heat equation (|2.12p on the rect- 
angular grid {(/tk, t£)} is given by the (n + l)-dimensional convolution 

f htT (hk,Ti) = (w+1)/2 r^- Yl K 2 (h(k-m),Te,Ti)u(hm,Ti) , 

meZ" 

where the integral if 2 given by (|2.14p cannot be taken analytically. Making the substitu- 
tion 

1 + e"« 

K2(x,t,ri) transforms to the integral over R 



t J e -(i/(l+C«)-Ti) 2 /(»0T 2 ) e -|x|2/(m2+4^/(l+C-«)) 

K 2 (*,t,Tl) = - j — - 4 ^ /(i - e _ f))n/2 cosh 2 (?/2) ^ 

— oo 

with exponentially decaying integrand. The quadrature rule for K2(hk, t£, ti) 



T ^ Wp (Vh? + + e-^)y/2 cos h 2 (e p /2) i = A 



e ^ 



fc2/(D+4^£/(h 2 (l+e-«f))) 



with certain weights and nodes £ p and the separated representation of the right-hand 
side of ([2TT2]) 

Q n 

u(hx,t) = j2r q l[u < f\x j ,t) + 0(e) 

9=1 3=1 



leads to the approximation 



f h>T (hk,r£) « 4ff („ +1 )/ 2 ^EE c J ( y 2) E^^llp- 



where we use the abbreviation 

e _(£ /( l +e « P )_ i )2 /7?0 

°^' p = (P/i 2 + W/(l + e-«p))™/ 2 ' 
and j is the product of one-dimensional convolutions 

n 

cg )i:P = n J] (fcm ij Ti)e-( fc '-^) a /(2'+^W a (i+e- e ''))) . 

j=l m,j&L 

Approximations which converge with higher order to the solution of (|2,12p can be 
obtained if the right-hand side u is approximated by 

UhA ^ t] = -jnw § u{hm > Ti)?}2S V7^) n^r^n (3 - 10) 

meZ" 
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with the basis functions r]2M defined by (|3.3|) . Then by using (|3.5p the heat potential of 
the quasi-interpolant (|3.10p is given by 

t 

/h n 

(Pt-\Uh !T (-,\))(x)d\ = , / —^ = y~] K s ,M(x-hm,t,Ti)u(hm,Ti), 

u m£Z" 

where we denote 

fc=o ' * j=i e=o ° o 



We have 



(-1)' gj? p -|x-te,|V(P* a -Hu>A) _ e-|x-^l 2 /(m 2+ 4,A) _ i/2 )( ( x .- hm .) 



- x-/im - (•/'/;-+ i;'A) _ ; r (-1/2) / V^j "- '"-j ^ 2 

4^ dxf (P/i 2 + 4z^A)^ ' I M 2 + 4i/A 



which by using ()2.14p leads to the representation with separable integrand 



KsM x,t, ri) = 2 <ft T » / .-<HMWM J! »< A ' ' 

fc=0 g j=l 

Here we denote by gu the function 

„„Q ~,\ —x 2 / (T>h 2 +4v\) ^ST~~^ {Vh 2 f .(-1/2)/ X 2 



4 Numerical results for the Newton potential 

We provide results of some experiments which show the accuracy and numerical conver- 
gence orders of the method. In the cube [—6, 6] n , n > 3, we compute the Newton potential 
of the densities 

Ul (x)=e-I x l 2 and u 2 (x) = (4|x| 2 - 2n) e"l x l 2 , 
which have the exact values 

£ n iti(x) = 4 | x | L w _2 ^(t; ~ 1; l x ! 2 ) and £ n u 2 (x) = -e" |x|2 . 

In Table ffl we compare the values of the exact and the approximate solution for C n u\ 
at some points (xi, 0, . . . , 0) G M n of the grid, where we have chosen the space dimensions 
n = 3,10,100,300. 

The approximations have been computed on a uniform grid with step size h = 0.05 
using the basis functions 

n 

7j»(x) = (TrVy n / 2 e~M 2 / v Yl L^ /2 \x 2 /V) with V = 3.5 . 

i=i 
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dimension 


X\ 


exact 


approximation 


relative error 






A A 

0.0 


a c a a a a a a a a a a aaa 

0.50000000000000 


a /Iaaaaaaaaaoocta 

0.49999999923850 


1.5230E-09 






1.0 


0.37341206640621 


0.37341206614375 


T AOO I 7T7' 1 A 

7.0287E-10 


n = 


= 3 


a A 

2.0 


A OA A £ O AO A A A/ 1 1 

0.22052034769061 


a ooacoao a ^ n n c\ a o 

0.22052034766043 


1.3685E-10 






o A 

3.0 


(J.1477l)lzz47l)99z 


A 1 yl^^AI 00^7A^ AO 

U.1477(Jlzz47(J4zo 


3.8549E-11 






A A 

4.0 


0. 1 1077836397370 


A 1 1 A^^OOrOA/'rCO 

0.11077836396658 


P A O /i ATP 1 1 

6.4242E-11 






5 n 


08869969954514 


08869969953834 


7 6764E-1 1 

1 .Ul U4:lJ J_ J_ 






0.0 


A A/ 1 A P" ACt A A ("in A A A A 

0.06250000000000 


A A/'O JAAAAAOOA^O 

0.06249999932963 


1 A^O/ 1 TT 1 AO 

1.0726E-08 






i a 

1.0 


A AAO A OAAO CT O 1 /iAO 

(J.Uz848zz353142o 


A A AO A O A AO E A /I E A A 

(J.(Jz848z2o5U459l) 


A A AAATT 1 AA 

9.4zl)9-h-U9 


n = 


10 


A A 

2.0 


A AAOOI ACT1 1 A1 O A O 

0.00331951101348 


A AAOOI AC1 AAA71 A 

0.00331951099712 


A AAOA17 AA 

4.9z8L)-h-U9 






O A 

3.0 


A AAAAAOTTAO ATO A 

0.00022377080789 


A AAAAA077A0A0 E 1 

0.00022377080851 


A TT A 1 "TT 1 A A 

2.7741_h-U9 






A A 

4.0 


A A A A AAAO 0£? A CT 1 T C 

(J.UUUUzz886U5175 


A AAAA A AO O A E 1 £3 A 

O.UU(J(J2288oU51o9 


A 1 ATTT 1 AA 

z\blz7-h-U9 






O.yJ 


u.uuuuuoooyjyyoi 


U.UUUUuOOOJc'Jc'OO 








A A 

0.0 


A A A r 1 AO A 1 AO 1 CO O 

0.00510204081633 


A A A E 1 A A A \ O O/"' 4 

0.00510204386664 


E ATO^fTT 1 AT 

5.9786E-07 






1 A 

1.0 


A AA1 A1 roOC1 AO 1 A 

U.Ul)1915zz51z312 


A AA1A1 CAAPAAA^A 

U.UU1915zzoz(Jz7z 


E COCAT7 AT 

5.6369E-07 


n = 


100 


A A 

2.0 


A AAA1A1 CTCOAOITA 

0.00010155802170 


A AAA 1 A 1 C COAOAAC 

0.00010155808096 


E O O /I TT7 1 AT 

5.8347E-07 






O A 

3.0 


A AAAAAATC71 A dO'7 

0.00000076714427 


A AAAAAAW71 /1EAO 

0.00000076714503 


A AAAA17 AT 

9.99z9hj-(J7 






A A 

4.0 


A AAAAAAAAAD A AO E 

0.00000000084085 


A AAAAAAAAAO A AO E 

0.00000000084085 


1 O O A 1 TT 1 f\f 

1.88Ulhj-(Jo 






o.u 


U.UUUUUUUUUUUUl^ 


n nnnnnnnnnnnni a 


O.D / UZJZj-UO 






0.0 


0.00167785234899 


0.00167786399030 


6.9382E-06 






1.0 


0.00062138979909 


0.00062139403983 


6.8246E-06 


n = 


300 


2.0 


0.00003157272440 


0.00003157294168 


6.8819E-06 






3.0 


0.00000022027431 


0.00000022027615 


8.3417E-06 






4.0 


0.00000000021134 


0.00000000021134 


8.4873E-06 






5.0 


0.00000000000003 


0.00000000000003 


2.6541E-05 



Table 1: Exact and approximated values of £ n ui(xi, 0, . . . , 0) and the relative error using £ 



In this case one has to determine 



4j Vl = l^ (1+t)* fe VX>(l + t)// (l+t)n/2 







for the integer vectors m = (mi, . . . , m n ) £ [—240, 240] n , which reduces by a quadrature 

W^^a+r^/ 2 ^^ (1 + 7*)* * ^(1 + 72) 
to the computation of one vector of length 481 with the components 

V ^- —L { - 1/2) — — -V m = -240,..., 240, (4.11) 

(i + r,) fc k Kv^i + n)) 1 ' ' ' v ; 

for any quadrature node and M = 4. Let us mention, that for all calculations, reported 
here and in the following Table O the same quadrature rule is used for computing the 
integrals Is(m/\^D). 

Table [2] shows that the cubature method is effective also for much higher space dimen- 
sions. We compare the exact values of C n U2 and the approximate values I^\u2 at the 
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n 




10 000 


100 000 


200 000 


X\ 


exact 


abs. error rel. error 


abs. error rel. error 


abs. error rel. error 


0.0 
1.0 
2.0 
3.0 
4.0 
5.0 


-1.0000E-00 
-3.6788E-01 
-1.8316E-02 
-1.2341E-04 
-1.1254E-05 
-1.3888E-11 


5.876E-05 5.88E-05 
2.160E-05 5.87E-05 
1.077E-06 5.88E-05 
7.345E-09 5.95E-05 
6.957E-12 6.18E-05 
9.304E-16 6.70E-05 


2.041E-03 2.04E-03 
7.509E-04 2.04E-03 
3.739E-05 2.04E-03 
2.522E-07 2.04E-03 
2.306E-10 2.05E-03 
2.863E-14 2.06E-03 


2.153E-03 2.15E-03 
7.920E-04 2.15E-03 
3.944E-05 2.15E-03 
2.659E-07 2.15E-03 
2.429E-10 2.16E-03 
3.008E-14 2.17E-03 



Table 2: Exact values of C n u%{xi, 0, . . . , 0), absolute and relative errors using £4 J 025 
same grid points (x\, 0, . . . , 0) £ R n for space dimensions n = 10 5 , 10 6 , 2 • 10 6 . The param- 

(n) 

eters of C\ / are T> = 3.5 and h = 0.025. The computing times on a Xeon processor with 
3.4 GHz are 0.42 seconds for space dimension n = 10 000, 5.12 seconds for n = 100 000, 
and 11.63 seconds for n = 200 000. 



M = 4 


ii 


3 


10 


500 


2 000 


30 000 




error rate 


error rate 


error rate 


error rate 


error rate 


5 
10 
20 
40 
80 


4.99E-05 
4.73E-07 6.72 
2.32E-09 7.67 
9.64E-12 7.91 
4.99E-14 7.60 


6.33E-04 
4.16E-06 7.25 
1.88E-08 7.79 
7.64E-11 7.94 
4.02E-13 7.57 


3.93E-02 
2.62E-04 7.22 
1.17E-06 7.81 
4.75E-09 7.94 
2.50E-11 7.57 


1.34E-01 
1.05E-03 6.99 
4.69E-06 7.81 
1.91E-08 7.94 
1.00E-10 7.57 


3.67E-01 
1.55E-02 4.57 
7.04E-05 7.78 
2.86E-07 7.94 
1.51E-09 7.57 


M = 3 


5 

10 
20 
40 
80 


1.45E-04 
5.05E-06 4.84 
9.76E-08 5.69 
1.61E-09 5.92 
2.55E-11 5.98 


4.11E-03 
9.35E-05 5.46 
1.62E-06 5.85 
2.60E-08 5.96 
4.09E-10 5.99 


1.98E-01 
6.23E-03 4.99 
1.08E-04 5.85 
1.73E-06 5.96 
2.72E-08 5.99 


3.51E-01 
2.44E-02 3.85 
4.34E-04 5.81 
6.95E-06 5.96 
1.09E-07 5.99 


3.68E-01 
2.37E-01 0.64 
6.46E-03 5.19 
1.04E-04 5.95 
1.64E-06 5.99 


M = 2 


5 
10 
20 
40 
80 


1.43E-03 
1.04E-04 3.78 
6.99E-06 3.90 
4.46E-07 3.97 
2.80E-08 3.99 


2.89E-02 
2.32E-03 3.64 
1.55E-04 3.91 
9.83E-06 3.98 
6.17E-07 3.99 


3.66E-01 
1.29E-01 1.51 
1.04E-02 3.63 
6.66E-04 3.96 
4.18E-05 3.99 


3.68E-01 
3.02E-01 0.28 
3.98E-02 2.92 
2.67E-03 3.90 
1.68E-04 3.99 


3.68E-01 
3.68E-01 0.00 
3.02E-01 0.28 
3.81E-02 2.99 
2.51E-03 3.92 


M = 1 


5 

10 
20 
40 
80 


3.73E-02 
9.29E-03 2.00 
2.31E-03 2.01 
5.75E-04 2.00 
1.44E-04 2.00 


1.93E-01 
6.56E-02 1.56 
1.79E-02 1.88 
4.56E-03 1.97 
1.15E-03 1.99 


3.68E-01 
3.68E-01 0.00 
3.51E-01 0.07 
1.99E-01 0.82 
6.50E-02 1.61 


3.68E-01 
3.68E-01 0.00 
3.68E-01 0.00 
3.52E-01 0.07 
1.99E-01 0.82 


3.68E-01 
3.68E-01 0.00 
3.68E-01 0.00 
3.68E-01 0.00 
3.68E-01 0.00 



Table 3: Absolute errors and approximation rates for £ n it2(l, 0, . . . , 0) using C 



In Table we report on the absolute errors and the approximation rates for the New- 
ton potential C n U2(l, 0, . . . , 0) = —0.3678794411714423 in the space dimensions n = 
3,10,500,2 000, and 30 000. The approximate values are computed by the cubature for- 
mulas C^M h defined by (|3.4p for M = 1, 2, 3, 4, having the approximation orders (|2.6p . with 
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T> = 5. We use uniform grids of size h = 0.1-2 1 - k , fc = 0,...,4, i.e., on the finest grid one 
has to compute vectors of length 1921 with the components (|4.1ip for m = —960, . . . , 960. 

The numerical results show that the high order cubature formulas give essentially bet- 
ter approximations of the Newton potential than second order formulas. For very high 

(n) 

dimensional cases the second order formula fails, whereas the 8-th order formula C\ I 
approximates with the predicted approximation rate. 



5 Quadrature rules for integral representations 

Here we describe the quadrature rule used to get separated representations of the volume 
applied to the basis functions. We give some numerical results on the quadrature error for 
second and fourth order formulas approximating the Newton potential and the potential 
of-A + a 2 . 



5.1 Quadratures 

It is well known that the classical trapezoidal rule is exponentially converging for certain 
classes of integrands, for example periodic functions or rapidly decaying functions on the 
real line. For any sufficiently smooth function, say of the Schwartz class <S(M), Poisson's 
summation formula yields that 

h E w) = E f^t) ■ 

k=—oo k=—oo 

Here / is the Fourier transform 

oo 

/(A)= J f(x)e~ 2mxX dx. 

— oo 

Thus, 



/ 



f(x)dx-h E f(kh) = J2f{- 



h 

k=~oo fc^O 



which indicates that by choosing special substitutions, which transform the integrand to a 
rapidly decaying function with rapidly decaying Fourier transform, the trapezoidal rule of 
step size h can provide very accurate approximations of the integral for a relatively small 
number of quadrature nodes {kh}. 

Additionally, if the integrand is analytic in a strip = {z £ C : | lmz\ < d} such that 

N(f,D d ) = \\f(- + id)\\ Lim + \\f(--id)\\ Ll{R) <oo, 

then results from Sine approximation can be applied to derive error estimates for the 
trapezoidal rule. It was shown in [HE], that for doubly exponentially decaying /, i.e., 

1/0*01 < Cexp (- ae i|x| ) for all x G M with constants a, b, C > , (5.12) 
the truncated rule with h = \og(2n aN /b) / (aN) satisfies 

/ f(x)dx-h E f( kh )\ < CN(f,D d )e- 27radN/los{27TaN/b) . 



-N 
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5.2 Newton potentials 



We make the substitutions 

i = e ? , £ = a(r + e r ) and r = b(u - e~ u ) , (5.13) 



with certain positive constants a, b, proposed in Waldvogel |15| . Then the integrals (|3.6|) 
are transformed to integrals over R with integrands /(u,x) decaying doubly exponentially 
in u. After the substitution we have 



oc 



ii(x) = ab 



(l + e-")(l + exp(6(u-e-")))0( M ) wa/(1+ * (u)) 
(l + 0(n))"/ 2 



with ^)(w) = exp(ab(u — e ") + aexp(6(u — e "))). Similarly 



/ (l + e-")(l + exp(b(^-e-")))0H fr 
J M (x) = a6 y 1 1 <?M («, ) du 

— oo -? = ^ 

with the function 

mM ..^iig 4-^Aiy.))) , 

For the practical application it is important that the quadrature rule 

N\ 

f(u,x)du*h ]T f( hk > x ) ( 5 - 14 ) 

k=N 

approximates /m( x ) with prescribed error uniformly for |x| < K with a minimal number 
of summands. This number strongly depends on K and on the parameters a, b in the 
transformation f)5. 13j) . For the computations in Section [H which cover a wide range of 
dimensions and different orders of the cubature formulas, we did not try to find optimal 
parameter sets. The numbers given in Tables [1] and [2] were obtained with the parameters 
a = b = 2 in (l5TT3|) and h = .02, N = -35, JVi = 80 in (f5TT4]h For the results of Table [3] 
we have chosen a = 6, b = 5 and h = .003, iVo = 39, N\ = 250. We report in the following 
on some test to determine optimal parameters for second and fourth order formulas in low 
dimensional cases (see also [TO]). 



5.2.1 Approximation to the integral Ii(x) 

We assume in (|5.13p a = b = 1. Figure [TJ illustrates the graph of the integrand function 
f(u, x), u G (—4,4), n = 3, for different values of |x| < 10 3 . A similar behavior holds for 
other dimensions n. Table H] presents the maximum step h§ and the minimum number of 
quadrature points required to achieve the relative error e, uniformly in |x| G [0, 10 3 ] for 
space dimensions n = 3, 4, 5, 6. 

It is possible to play with different parameters a and b in order to diminish the number 
of summands in the quadrature formula. Consider e.g. the case a = 6 and 6 = 5. Figure [5] 
shows the graph of f(u, |x|), u G (0,0.85) for different values of |x|. The numerical results 
for this quadrature are given in Table [5j 
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1 



-4 -2 2 4 -4 -2 2 4 -4 -2 2 4 

Figure 1: The plot of the integrand function f(u,x) (a = b = 1) in Ii(x) for |x| : 
0, 100, 1000 (from the left to the right) in the interval u £ (-4,4). 





n = 3 


n = 4 


n = 5 


n = 6 


rel. error 


h 


nodes 


h 


nodes 


h 


nodes 


h 


nodes 


1E-05 


0.072 


82 


0.072 


77 


0.065 


83 


0.059 


90 


1E-07 


0.055 


116 


0.051 


121 


0.060 


96 


0.044 


130 


1E-09 


0.043 


161 


0.040 


164 


0.037 


169 


0.035 


178 


IE- 11 


0.036 


205 


0.033 


206 


0.033 


200 


0.029 


220 



Table 4: The approximation of Ii(x) for |x| < 10 3 , with a = b = 1 in (|5.13[) 





n — 3 


n = 4 


n — 5 


n = 6 


rel. error 


ho 


nodes 


h 


nodes 


h 


nodes 


h 


nodes 


1E-05 


0.0077 


61 


0.0070 


83 


0.0069 


57 


0.0058 


70 


1E-07 


0.0055 


111 


0.0049 


96 


0.0046 


112 


0.0042 


117 


1E-09 


0.0042 


170 


0.0037 


169 


0.0034 


179 


0.0037 


158 


IE- 11 


0.0034 


247 


0.0033 


200 


0.0031 


221 


0.0028 


242 



Table 5: The approximation of Ji(x) for |x| < 10 3 , with a = 6, b — 5 in (|5.13|) 



5.2.2 Approximation to the integral l2(x) 
Next, we discuss the computation of the integral 

00 n 2 

hU) - f T\c~ x ' /(t+1) ( 3 + 2t - Xj ) - 

m*j -y ne ^ j (1+t)n/2 

o J =1 

using the variable transformation (|5.13p and the trapezoidal rule (|5,14p . for n = 3 and 
n = 4. The graphs of the integrands f(u,x) are very similar to the case of ii(x). In the 
numerical results, for the sake of simplicity, we assumed x = (x, x, x), with |x| < 10 3 . The 
maximum step ho and the minimum number of quadrature points required to achieve the 
relative error e uniformly in |x| G [0, 10 3 ] are presented in Table [6] for the cases a = b = 1 
and a = 6, b = 5. 



16 



30 



20 



10 




0.4 
0.3 
0.2 
0.1 



0.3 0.4 0.5 0.6 0.7 0.8 




0.05 
0.04 
0.03 
0.02 
0.01 



0.3 0.4 0.5 0.6 0.7 OA 



0.3 0.4 0.5 0.6 0.7 0.£ 



Figure 2: The plot of the integrand function /(it, x) (a = 6, b = 5) in ii(x) for |x| 
0, 100, 1000 (from the left to the right) in the interval u G (0,0.85). 





n = 3 


n = 4 


n = 3 


n = 4 


rel. error 


h 


nodes 


h 


nodes 


ho 


nodes 


h 


nodes 


1E-05 


0.072 


82 


0.072 


77 


0.0077 


63 


0.0074 


57 


1E-07 


0.055 


118 


0.051 


121 


0.0052 


114 


0.0046 


120 


1E-09 


0.043 


163 


0.040 


163 


0.0042 


175 


0.0037 


175 


IE- 11 


0.036 


204 


0.033 


206 


0.0034 


234 


0.0033 


222 



Table 6: The approximation of /2(x) for |x| < 10 3 , with a = b = 1 (left) and a = 6, b = 5 (right) 



5.3 Potentials of advection-diffusion operators 

We consider the special case b = 0, c = a 2 > 0, which leads to the one-dimensional 
representation of the volume potential applied to n r\iM 

JL ,^£=1 «-*?/(!+*) , r 2 ^ dt 

M > 1 



(cf. (|3.8p . (|3.9p ). To get doubly exponential integrands we make the substitution 

t = <p(u) with 4>{v) = exp(b(u — exp(— u)) , b > , (5.15) 
and apply the trapezoidal rule to 

y + 11 2-; (1 

Figures [3] and [3] illustrate the graph of the integrand /(it, x) of ifi(x), n = 3, for different 
values of |x| and a 2 = 0.01 and a 2 = 4. 

In Table [7] we present the maximum step /to and the minimum number of quadrature 
points in formula (|5. 14|) for the integral -f^i(x) for n = 3, in the cases a 2 = 0.01,0.1, 1,4. 
We have chosen b = 1 in the substitution (j5.15|) . 
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0.8 
0.6 
0.4 

V 



8 10" 6 
6-10" 6 
4-10" 6 
2-10" 6 



2.5 10" 45 
2-10" 45 

1.5-10" 45 
1-10" 45 
5 10" 46 



Figure 3: The plot of the integrand function f(u, x) in -Ki(x), a = 0.01,5 = 1 for 
|x| = 0, 100, 1000 (from the left to the right) in the interval u £ (-3, 3). 




6-10" 



4-10" 



2 10" 44 



2-10" 174 
1.5-10" 174 
1-10" 174 
5-10" 175 



Figure 4: The plot of the integrand function f(u, x) in K"i(x), a 2 
0,50,200 (from the left to the right) in the interval u £ (-2,2). 



4,6 



1 for Ixl = 



a' 


0.01 


0.1 


1 


4 


rel. error 


ho 


nodes 


ho 


nodes 


^0 


nodes 


ho 


nodes 


1E-05 


0.58 


20 


0.58 


17 


0.48 


15 


0.44 


13 


1E-07 


0.47 


25 


0.42 


16 


0.38 


20 


0.36 


17 


1E-09 


0.39 


32 


0.40 


25 


0.37 


22 


0.31 


21 


1E-11 


0.30 


43 


0.29 


36 


0.29 


28 


0.27 


25 


1E-13 


0.26 


50 


0.25 


43 


0.25 


34 


0.25 


29 



Table 7: The approximation of ^i(x) with different values of a 2 for |x| < 10 3 
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